function [rej_rate,p_adj] = rej_rate_CCE_mean(SigmaHat,membership,alpha_sig,alt)

rng(7)
n = size(SigmaHat,1);
clusters = unique(membership);
G = numel(clusters);

nSim = 100000;  % number of simulation replications

% only generate group mean
L = repmat(membership',G,1)==repmat(clusters,1,n);
N_g_vec = sum(L,2);
L_weighted = L./repmat(N_g_vec,1,n);
V = alt+(randn(nSim,G)*chol(L_weighted*SigmaHat*L_weighted'))';
mu_hat_vec = (sum(L,2)/n)'*V; % mean estimates
LU_hat = repmat(N_g_vec,1,nSim).*V-L*ones(n,1)*mu_hat_vec; % L*U_hat
se_vec = sqrt(sum((LU_hat).^2)/n^2); % summing up related correlations
t_vec = mu_hat_vec./se_vec;
test_vec = abs(t_vec)>tinv(1-alpha_sig/2,G-1);
rej_rate = mean(test_vec);

% adjusted critical values and corresponding p-value threshold such that 
% rejection rate is alpha_sig
cv_adj = quantile(abs(t_vec),1-alpha_sig); 
p_adj = 2*(1-tcdf(cv_adj,G-1));
